module sw_case5_mod

  use const_mod

  implicit none

  private

  public sw_case5

  real(8), parameter :: radius  = 6371.220d3
  real(8), parameter :: gravity = 9.80616d0
  real(8), parameter :: omega   = 7.292d-5
  real(8), parameter :: u0      = 20.0d0
  real(8), parameter :: gh0     = 5960.0d0 * gravity
  real(8), parameter :: ghs0    = 2000.0d0 * gravity
  real(8), parameter :: R       = pi / 9.0d0
  real(8), parameter :: lon0    = pi * 1.5d0
  real(8), parameter :: lat0    = pi / 6.0d0
  real(8), parameter :: alpha   = 0.0d0
  real(8), parameter :: cos_alpha = cos(alpha)
  real(8), parameter :: sin_alpha = sin(alpha)

contains

  subroutine sw_case5(lon, lat, u, v, gh, ghs, g, f, radi)

    real(8), intent(in ) :: lon
    real(8), intent(in ) :: lat
    real(8), intent(out), optional :: u
    real(8), intent(out), optional :: v
    real(8), intent(out), optional :: gh
    real(8), intent(out), optional :: ghs
    real(8), intent(out), optional :: g
    real(8), intent(out), optional :: f
    real(8), intent(out), optional :: radi

    real(8) cos_lon, sin_lon, cos_lat, sin_lat, dlon, d
    
    if(present(radi))radi = radius
    
    cos_lon = cos(lon)
    sin_lon = sin(lon)
    cos_lat = cos(lat)
    sin_lat = sin(lat)
    dlon    = abs(lon - lon0)
    dlon    = min(dlon, 2 * pi - dlon)
    d       = min(R, sqrt(dlon**2 + (lat - lat0)**2))

    if (present(u  )) u   = u0 * (cos_lat * cos_alpha + cos_lon * sin_lat * sin_alpha)
    if (present(v  )) v   = -u0 * sin_lon * sin_alpha
    if (present(gh )) gh  = gh0 - (radius * omega * u0 + u0**2 * 0.5d0) * (sin_lat * cos_alpha - cos_lon * cos_lat * sin_alpha)**2
    if (present(ghs)) ghs = ghs0 * (1.0d0 - d / R)
    if (present(g  )) g   = gravity
    if (present(f  )) f   = 2 * omega * (-cos_lon * cos_lat * sin_alpha + sin_lat * cos_alpha)

  end subroutine sw_case5

end module sw_case5_mod
